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The  work  to  date  has  dealt  with  the  effect  of  seismometer  placement  on  array 
sidelobes,  the  performance  to  be  expected  with  arrays  using  one  bit  signals, 
from  the  seismometer,  the  coherence  of  signals,  noise,  and  coda  across  LASA, 

P-ooda  correlation  as  a  discriminant,  and  the  frequency  dependence  of  attenuation 
in  the  earth's  upper  mantle. 

A  randomized  array,  as  Opposed  to  ono  with  clustered  elements  has  show 
isdolobeo  3ome  14  dR  lower  than  those  of  the  current  LASA  design.  The  one  bit 
array  (DIMUS)  has  shown  a  processing  gain  of  2  dR  less  than  that  for  an  analog 
array,  no  predicted  by  theory.  The  same  array  with  the  addition  of  a  singlo 
analog  channel  offers  good  posaibll itios  of  proving  as  useful  in  identification 
a3  a  conventional  array. 

Main  P,  and  noiee  correlations  observed  across  LASA  Montanna  from  LONGSfOT 
showed  signal  correlation  to  bo  high,  0,97,  and  noise  quite  low  for  spaclngo  equal 
to  thoso  of  the  cluotoro.  P-coda  correlation  decreased  to  a  valuo  noar  zero  for 
soiomomotoi  spacings  on  the  order  of  150  km. 

Study  of  F -coda  correlation  coefficients  for  widoly  spaced  seismometers 
fell  into  two  distinct  distributions  with  an  equivalent  f/n  of  1A  DB. 

A  method  wa3  developed  for  estimating  the  Q  of  the  upper  mantle  using  the  P 
end  pP  phases  of  deep  earthquakes.  First  test  of  the  method  indicate  an  average 
Q  of  approximately  200  for  the  upper  500  Km  of  tho  mantle  in  the  Fiji  Islands 
Region, 


DD  Mil,  1473 


Security  Classification 


GENERAL  ATRONICS  CORPORATION 


5  July  1966 


Report  1531-2033-7 
SYNOPSIS 

This  report  covers  the  work  performed  during  the  first 
quarter  of  Contract  AP19 (626) -5961,  Large  Aperture  Seismic  Arrays. 

•  The  work  to  date  has  dealt  ,/ith  the  effect  oi  seismometer  place¬ 
ment  on  array  sidelobes,  the  performance  to  be  expected  with 
arrays  using  one  bit  sigrals  from  the  seismometers,  the  coherence 
of  signals,  noise,  and  coda  across  LASA,  P-coda  correlation  as 
a  discriminant.  and  the  frequency  dependence  of  attenuation  in 
the  earth's  ujpsr  mantle. 

A  randomized  array,  as  opposed  to  one  with  clustered  elements, 
has  shown  sidelobes  some  14  dB  lower  than  those  of  the  current 
LASA  design.  The  one  bit  array  (DIMUS)  has  shown  a  processing 
gain  of  2  dB  less  than  that  for  an  analog  array,  as  predicted 
by  theory.  The  same  array  with  the  addition  of  a  single  analog 
channel  offers  good  possibilities  of  proving  as  useful  in  identi¬ 
fication  as  a  conventional  array. 

Main  P,  and  noise  correlations  observed  across  LASA  Montanna 
from  LONGSHOT  showed  signal  correlation  to  be  high,  0.97,  and 
noise  quite  low  for  spacings  equal  to  those  of  the  clusters. 

P-coda  correlation  decreased  to  a  value  near  zero  for  seismo¬ 
meter  spacings  on  the  order  of  150  km. 

Study  of  P-coda  correlation  coefficients  for  widely  spaced 
seismometers  fell  into  two  distinct  distributions  with  an 
equivalent  S/N  of  16  dB. 

A  method  was  developed  for  estimating  the  Q  of  the  upper 
mantle  using  the  P  and  pP  phases  of  deep  earthquakes.  Pirst 
test  of  the  method  indicate  an  average  Q  of  approximately  200 
for  the  upper  500  km  of  the  mantle  in  the  Fiji  Islands  Region. 

This  research  was  sponsored  by  the  Advanced  Research  Pro¬ 
jects  Agency  and  was  monitored  by  ESD  (ESL)  under  Contract 
Number  AP19 (628) -5961. 
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Section  I  Introduction 

This  is  the  first  quarterly  technical  report  on  Contract 
AE19  (626)-5961,  Large  Aperture  Seismic  Arrays.  The  subject 
matter  of  this  program  represents  extensions  to  earlier  work 
in  connection  with  the  systems  aspects  of  the  detection  and 
identification  of  underground  explosions  performed  under  Con¬ 
tracts  SD  251  and  SD  283.  Work  on  the  current  program  deals 
with  five  general  topless 

1)  Array  Theory 

2)  Enhancement  of  pP  in  coda 

3)  Application  of  linear  discriminant  analysis  to 
short  range  discriminants 

4)  High  frequency  content  of  seismograms 

5)  Utility  of  signal  processing  against  tests  con¬ 
ducted  so  as  to  resemble  or  be  masked  by  natural 
events . 

Included  in  the  first  topic  are  not  only  theoretical  determine 
tion  of  beam  properties,  but  also  measurements  of  signal  and 
noise  properties  which  affect  beam  formation. 

Work  during  this  first  quarter  was  concentrated  on  the 
first  and  fourth  topics,  with  most  emphasis  placed  on  the 
first.  Here  the  effects  of  the  distributions  of  sensors  on 
array  performance  foi  teleseismic  signals  was  considered  - 
particularly  with  respect  to  the  possibilities  of  control  of 
array  sidelobes .  Also  studied  were  the  correlation  of  signals 
noise,  and  P-coda  with  distance  in  so  far  as  they  would  affect 
performance  of  arrays.  An  offshoot  of  this  work  was  the  dis¬ 
covery  that  the  zero  offset  coda  correlation  coefficient  is  a 
potentially  good  discriminant  for  separation  of  shallow  earth¬ 
quakes  and  underground  explosions.  Finally,  investigation  was 
made  of  the  effects  of  array  performance  when  only  the  zero 
crossings  of  the  individual  instrument  outputs  is  retained. 
Here,  potentially  large  savings  in  data  storage  and  trans¬ 
mission  equipment  may  be  realizable.  Simpler  signal  process¬ 
ing  equipment  may  also  be  possible  where  these  simpler  repre¬ 
sentations  of  the  signals  can  be  used.  Moreover,  it  may  be 
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possible  to  avoid  clustering  of  seismometers  -  thereby  improv¬ 
ing  the  sidelobe  characterisitcs  of  arrays,  if  it  is  possible 
to  use  hard  limited  signals  from  the  individual  seismometers. 
To  date,  this  investigation  shows  promise. 

Work  on  the  high  frequency  aspect  of  seismic  signals  has 
been  limited  to  development  of  a  method  of  determining  attenu¬ 
ation  in  the  upper  mantle.  Results  of  further  work  with  this 
method  should  enable  choice  of  those  sites  particularly  suited 
to  reception  of  high  frequency  signals,  should  they  exist. 
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Section  II  Beam  Theory 

Section  II-A  DERIVATION  OF  THE  DIMENSIONS  OE  THE  TELES EISMIC 

BEAM  OE  A  LARGE  APERTURE  ARRAY 


In  earlier  reports  a  presentation  of  teleseismic  beam 
1  2 

theory  was  given.  *  Among  the  subjects  discussed  were  the  cross- 
sectional  properties  of  a  teleseismic  beam,  the  depth  of  field  of 
the  array  and  a  comparison  of  array  properties  for  several  seismic 
phases.  Given  also  were  a  radial  pattern  showing  side  lobes  from 
about  40°  to  100°  and  the  effects  upon  the  sidelobe  structure  of 
a  transient  waveform. 

Omitted  at  that  time  were  the  derivations  of  the  three  dimen¬ 
sions  of  a  teleseismic  beam.  These  derivations  are  the  subject 
of  this  section  of  the  report. 

Eigure  132  shows  an  array  located  at  the  pole  of  a 
spherical  earth.  Two  perpendicular  meridians  are  drawn. 

The  beam  is  focussed  at  point  P  which  is  located  on  one 
of  the  meridians  at  latitude  tt/2  -  A  .  Thi3  means  that 
the  elements  in  the  array  have  time  delays  built  into 
them  which  exactly  correct  for  the  differential  travel 
times  from  P  to  the  array  elements.  The  meridian  contain¬ 
ing  P  is  called  the  radial  meridian.  The  other  meridian 
is  called  the  transverse  meridian. 

Shown  also  is  a  distorted  ellipse  containing  P.  This 
ellipse  is  an  array  gain  contour.  The  radial  extent  of 
the  beam  is  twice  PR,  the  size  of  the  ellipse  on  the 
radial  meridian.  The  transverse  extent  of  the  beam  is 
twice  PT,  the  size  of  the  ellipse  on  the  latitude  line 
containing  P. 

This  constant  gain  contour  is  the  intersection  of  a 
distorted  constant  gain  ellipsoid  with  the  surface.  The 
"vertical"  axis  is  along  the  seismic  path.  The  depth 
of  field  of  the  array  is  twice  PZ,  the  extent  of  the 
beam  from  the  focal  point  to  the  constant  gain  contour 
along  the  seismic  ray  path. 
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For  a  surface  focus  half  of  the  distance  covered  by 
the  depth  of  field  is  outside  the  earth  and  therefore 
not  pertinent.  Thus  the  depth  of  field  could  be  very 
sensibly  defined  as  the  distance  PZ ,  rather  than  twice 
it.  However,  for  deep  focal  points  the  alternate  argu¬ 
ment  would  prevail.  For  generality,  therefore,  the 
depth  of  field  is  taken  as  twice  PZ.* 


Parameters 


T 


Y 

Y 
n 


t 

r 


9 

r 

Ax 


Ay 


Ah 

T 


A 

A 


h 

h 

T 


o 

A 


=  dominant  period  of  the  seismic  signal  (sec) 

=  transverse  extent  of  the  array  (radians) 

=  radial  extent  of  the  array  (radians) 

=  aperture  efficiency 

=  transverse  extent  of  beam' within  3  db  contour  (radians) 

=  radial  extent  of  beam  within  3  db  contour  (radians) 

=  transverse  extent  of  the  beam  within  3  db  contour  (km) 
li 

2x10  Qfc 

TT 

=  radial  extent  of  the  beam  within  3  db  contour  (km) 

2xl0^9t 

7T 

=  3  db  depth  of  field  (km) 

=  travel  time  (sec) 

circle  distance  variable  (radians) 

of  A  from  array  center  to  beam  focal  point  (radians) 
variable  (km) 
of  focal  point  (km) 

( sec/rad ) 


=  great 
=  value 
=  depth 
=  depth 
3T 


*The  above  four  paragraphs  are  quoted  from  Reference  1,  pp .  l^O-l^l. 


-4- 


GENERAL  ATRONICS  CORPORATION 


T 

T 


A  A 

hA 


2 

d  T 


9  A 
2 

9^T 


9  A  9h 


( sec/rad^ ) 

(sec/rad  km) 


A  ,h 
o  *  o 


Summary 
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(rad) 
( rad) 
(km) 


(9) 

(13) 

(15) 


Dependence  of  Beam  Dimensions  Upon  Array  Dimensions 

The  width  of  the  main  lobe  of  the  radiation  pattern  in 
any  plane  containing  the  beam  axis  is  dictated  by  the  length  of 
the  projection  of  the  array  onto  the  plane.  By  analogy  with  the 
case  of  the  homogeneous  medium,  imagine  that  the  beam  path  shown 
in  Figure  1  is  straightened  so  that  the  beam  cross-section  to  be 
derived  (containing  P,  R  and  T)  falls  in  the  plane  of  the  equator 
Then  from  the  previous  statement  it  may  be  seen  that  the  trans¬ 
verse  extent  in  the  beam  2PT  is  dictated  by  the  transverse  size 
of  the  array  20A,  which  is  the  extent  of  the  array  along 
the  transverse  meridian.  Similarly,  the  radial  size  of  the  array 
20B,  dictates  the  radial  dimension  of  the  beam  2PR.  Less  obvious 
is  the  fact  that  the  radial  extent  of  the  array  also  dictates  the 
depth  of  field,  2PZ. 
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Transverse  Beamwidth 

Figure  2  shows  the  pertinent  geometry.  The  angular  dis¬ 
tance  A  between  point  A,  an  angular  distance  y  from  array  center 
in  the  transverse  direction,  and  point  Q,  an  angular  distance  9 
from  beam  axis  along  the  transverse  dimension  of  the  beam  cross- 
section,  is  given  by 

cosA  =  cosacos  6cos  l  +  sinotsing 

=  sinysinA^sinO  +  cosycosA.  (1) 

a  and  g  are  the  latitudes  of  array  point  A  and  beam  point  Q. 

a  is  the  difference  in  longitude.  Aq  is  the  distance  from  array 

center  0  to  beam  center  P. 

Subtraction  of  cosA  from  both  sides  gives 

o 

cosA-cosA  =  sinysinA  sinO  -  cosA  (1-cosy)  (2) 

o  '  o  o  ' 

The  difference  of  two  cosines  where  the  angular  difference  is 

small  is  well  approximated  by 

cosA  — cosAq  “  sinAQ • 6 A  (3) 

from  which 

6  A  -  sinysinO  -  ctnAQ ( 1-cosy )  (  4  ) 

(4)  represents  the  differential  distance  between  array  center 
to  beam  center  (Aq)  and  the  distance  from  point  A  in  the  array 
to  point  Q  in  the  beam  (A).  The  corresponding  differential 
arrival  time  is 

6T  =  6A  (5) 

where  3T/3A  is  the  first  partial  derivative  of  travel  time  with 
respect  to  distance.  If  the  outputs  of  the  seismometers  at  all 
points  of  the  array  were  added  simultaneously,  signals  arriving 
from  Q  would  be  separated  in  time  by  6T. 

To  form  a  beam  focussed  at  P  it  is  necessary  to  add  (or 
subtract)  a  time  delay  at  each  array  element  equal  and  opposite 
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to  (5),  so  that  the  seismic  energy  from  P  would  appear  simultaneously 

at  the  output  of  each  array  element.  The  correction  is  exactly  the 

second  term  in  (M),  since  the  first  term  is  zero  for  9  =  0.  This 

correction  is  a  function  only  of  the  focal  distance  A  and  the  trans- 

o 

verse  distance  A  of  the  seismometer  from  the  array  center.  It  is 
thus  a  permissable  correction  which  may  be  built  into  a  beamforming 
network.  Assuming  that  this  static  correction  is  made,  the  corrected 
differential  distance  becomes 

6A  =  sinysinO  (6) 


Combining  (5)  and  (6),  and  replacing  the  sine  by  the  angle,  yields 
the  very  simple  first  order  residual  differential  time  delay  after 
correction 


6T 


y  ,9<l/2 


(7) 


o 

Let  t  be  the  dominant  period  of  the  seismic  signal. 
Define  the  transverse  beamwidth  9^  as  the  range  in  9  over  which 
6T  varies  from  +t/4' to  -t/h.  This  definition  is  consistent  with 
the  usual  beamwidth  definition  of  a  uniformly  illuminated  array. 
The  first  step  is  (8)  in  which  the  notation- for  the  partial  deri¬ 
vative  has  been  changed  to  that  given  in  the  glossary. 


(8) 


(8)  gives  the  transverse  extent  of  the  beam  in  angular 
longitude  units  along  latitude  0 .  To  convert  this  quantity  to 
the  distance  or  angular  measure  on  the  surface  of  the  earth,  (8) 
must  be  multiplied  by  sinAQ.  Two  other  alterations  to  (8)  are 
worth  making.  First,  the  quantity  2y  in  the  denominator  may  be 
identified  with  the  transverse  extent  of  the  array,  for  which 
the  simpler  quantity  y,  is  substituted.  Secondly,  the  beamwidth 
will  normally  be  larger  than  (8)  due  either  to  amplitude  weight¬ 
ing  or  spatial  taper  in  the  aperture,  normally  dictated  by  side- 
lobe  tolerance  or  by  economics.  This  factor  is  usually  identified 
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with  aperture  efficiency  n,  a  quantity  typically  on  the  order  of 
2/3  to  3/^-  The  final  expression  is 


t  sinA 


VTA 


(rad) 


(9) 


All  angles  are  in  radians  and  all  time  quantities  are  in  seconds. 


Radial  Beamwidth 

Figure  3  shows  the  relevant  geometry .  A  point  in  the 

array  is  shown  a  distance  y  from  array  center  along  the  radial 

meridian.  The  axis  of  the  seismic  beam  is  shown  as  a  solid  line. 

The  focal  point  of  the  beam  is  a  distance  h  below  the  surface 

o 

of  the  earth.  A  point  in  the  radial  pattern  of  the  beam  a  distance 
9  from  the  focal  point  is  shown.  The  great  circle  distances  to 
the  array  center  from  this  point  and  from  the  focal  point  are  aq-9 
and  Aq,  respectively.  Let  T(A0,hQ)  be  the  travel  time  from  array 
center  to  the  focal  point  Ao,hQ  of  the  beam.  Let  T(Ac-Y,hQ)  be  the 
travel  time  from  the  point  in  the  array,  which  is  a  radial  dis¬ 
tance  y  from  the  array  center,  to  the  focal  point.  The  difference 
is  exactly  the  static  correction  for  the  seismometer  at  distance 
Y  when  the  array  is  focussed  at  AQ,ho.  This  quantity  is  (10). 

ST(Y3A0,ho)  =  T<Vho)  - 


9T 
9  A 


(10) 


A  ,  h 
o  o 


For  a  point  in  the  beam  a  distance  9  from  beam  axis 
along  the  radial  direction,  the  differential  travel  time  is 


6T(Y,A0-9,h0)  =  T(A0-9,hQ)  -  T(A0-Y-9,hQ)  (11) 

A  large  portion  of  the  differential  travel  time,  however,  is  re¬ 
moved  by  delay  correction  (10).  Thus  the  residual  differential 
travel  time  is  (11)-(10): 
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«T(y,Ve,h0)  =  TU0-@,h0)  -  T(40-Y-8,h0)  -  y  § 


=  Y 


3T 
3  A 


-  Y 


A  ,h 
o 5  o 


3T 
3  A 


A  ,h 

o  ’  o 


*  Y© 


p 

3^T 


3  A 


(12) 


A  ,  h 
o  ’  o 


Again,  identifying  the  radial  beamwidth  9^  with  that  range  in  beam 
angle  9  for  which  the  uncorrected  differential  travel  time  varies 
by  t/2,  identifying  2y  as  the  radial  extent  of  the  array  Yrj  and 
introducing  the  aperture  efficiency  factor  n,  gives  the  radial 
beamwidth . 


9 


V^i 


( rad) 


(13) 


Depth  of  Field 

The  geometry  of  Figure  3  suffices  for  this  derivation. 

The  residual  differential  travel  time  after  delay  correction  (10)  is 


6T(Y,Ao,h)  =  T(AQ,h)  -  T(Ao-Y,h)  -  Y 


A  ,h 
o 5  o 


3T 
3  A 


-  Y 


A  ,h 
o  ’ 


3T 
3  A 


A  h 
o  ’  o 


2 

.  y  3  ~T 
Y^h  ho^  3  A  3h 


(14) 


A  ,h 
o 5  o 


Folowing  the  earlier  procedure,  the  variation  in  depth  around 
the  focal  depth  hQ  over  which  the  residual  differential  travel 
time  varies  by  1/2  the  dominant  period  of  the  seismic  event  is 
taken  as  the  depth  of  field.  Thus,  the  depth  of  field  due  to  a 
radial  array  of  length  Yr  having  aperture  efficiency  t  is 
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Ah  »  -  T  ■  (km)  (15) 

Yr  1  hA 

The  derivation  of  (15)  involved  only  the  radial  extent 
of  the  array.  Although  it  is  not  immediately  obvious  it  is  only 
the  radial  length  which  has  significant  Influence  upon  the  depth 
of  field.  Hence  it  is  adequate  to  calculate  depth  of  field  in 
terms  of  the  radial  extent  only. 

To  demonstrate  this  fact  attention  is  called  to  (1*0. 

The  quantity  y  in  that  expression  represents  the  differential 
distance  between  a  point  in  the  beam  and  the  two  array  points 
along  the  radial  direction  in  the  array.  The  generalized  statement 
equivalent  to  ( 1 4 )  is  as  follows:  given  any  two  points  in  the 
array  y  in  (1*4)  is  replaced  by  the  distance  between  a  point  in 
the  beam  and  the  two  array  points.  In  particular,  if  the  array 
points  are  on  the  transverse  meridian,  it  is  necessary  to  return 
to  ( *4 )  ,  evaluated  for  9  =  0.  Thus  the  differential  distance 

6  A  =  ctnAQ ( 1-cosy )  (16) 

Since  y<<l 

6A  *  (y2/2)  ctnAQ  (17) 

the  transverse  equivalent  of  (1*4)  becomes 

2 

6T9(y,Ao,h)  =  2~  ctnAo(h-ho)ThA  (l8) 


Identifying  2y  =  y^  =  transverse  extent  of  the  array,  and  following 
the  same  procedure  as  before,  the  depth  of  field  due  to  the  trans¬ 
verse  extent  is 


AhQ  = 


*4x  tanA( 
Y02^ThA 


(19) 


If  the  array  is  symmetrical,  such  that  yy  =  y  ,  the  depth  of 
field  due  to  the  transverse  extent  is  greater  than  that  due  to 
the  radial  extent  by  the  ratio  of  (19)  to  (15),  or 


Ah 


9 


*4tanA 


Ah 


>>  1 


(20) 
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Thus  the  depth  of  field  of  the  array  is  dictated  by  its  radial 
size . 

As  an  example  consider  a  200  km  symmetrical  aperture 
focussed  at  Aq=  70°,  an  aperture  efficiency  of  3/4  and  a  domin¬ 
ant  seismic  period  of  one  second.  The  depth  of  field  due  to 
the  radial  dimension  is  3060  km.  If  the  array  were  collapsed 
to  a  line  array  in  the  transverse  direction,  the  depth  of 
field  would  increase  to  more  than  10^  km.  Even  for  very  large 
arrays  and  shorter  ranges  the  focal  point  remains  in  the  far 
field.  Thus  the  transverse  extent  of  the  array  does  not 
measurably  contribute  to  its  depth  of  field. 

The  three  beam  dimensions,  (9),  (13) ,  and  (15),  are 
plotted  in  Figure  4  for  the  P-wave.  This  figure  is  taken  from 
Reference  1,  Figure  2.  Travel  time  data  are  from  Jeffreys  and 
Bullen  .  In  the  figure  A  is  the  independent  variable.  The  de¬ 
pendent  variable  contains  the  beam  size,  the  effective  array 
size  Yh  and  the  dominant  period  t .  For  convenience  the  vertical 
scale  is  in  terms  of  km  degrees  rather  than  square  radians. 

Its  units  are  km  deg/sec.  To  obtain  beam  size  in  km  multiply 
the  ordinate  by  the  dominant  period  in  seconds  and  divide  by 
the  effective  array  size  in  degrees.  For  a  one  degree  array 
and  one  second  instruments,  the  vertical  scale  gives  the  beam 
dimensions  directly  in  km.  The  depth  of  field  curve  pertains  to 

a  surface  focus,  h  =0.  Because  of  the  manner  in  which  the 

o 

depth  of  field  is  defined  the  3  dB  contour  extends  below  the 
surface  only  to  one-half  the  depth  of  field.  Corresponding 
curves  for  any  other  seismic  phase  may  be  made  by  using  the 
travel  time  derivatives  for  that  phase. 
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Section  II-B  ARRAY  PATTERN  COMPUTATIONS 
General 

Previously  reported  studies  of  array  patterns  for  large 
seismic  arrays  have  been  primarily  concerned  with  the  perfor¬ 
mance  of  such  arrays  against  noise.  This  work  has  been  prin¬ 
cipally  performed  in  "k"  space,  and  as  a  result,  does  not  com¬ 
pletely  describe  the  array  performance  with  respect  to  tele- 
seismic  signals.  Those  areas  where  additional  information  can 
be  usefully  introduced  have  to  do  with  the  spreading  loss  and 
attenuation  of  seismic  propagation,  and  also  with  the  transient 
nature  of  seismic  signals.  Inclusion  of  these  effects  in  array 
pattern  calculations  give  rise  to  better  prediction-  of  the 
array  sidelobes  structure.  While  sidelobes  may  previously  have 
been  of  little  interest,  the  increased  signal-to-noise  gain 
achievable  with  arrays  such  as  LASA  may  convert  the  teleseismic 
noise  source  from  an  homogeneous  one  to  several  moving  and  time 
varying  "hot  spots"  representing  small  events  from  seismically 
active  regions  of  the  earth.  Hence,  it  becomes  important  to 
achieve  good  control  of  the  sidelobe  structure  of  these  large 
arrays . 

Results  to  date  do  not  take  station  correction  factors  into 
account,  particularly,  the  dependence  of  these  factors  on  range. 
It  is  believed,  however,  that  it  is  instructive  to  first  study 
the  case  where  such  corrections  are  assumed  to  be  constant. 

Later  work  can  extend  the  analysis  to  observed  variations. 

A  computer  program  taking  these  phenomena  into  account, 
was  developed  to  find  the  radial  and  transverse  beam  patterns 
of  seismic  arrays  for  P  wave  signals.  The  mathematical  model 
used  is  the  following: 

a)  The  array  seismometers  are  identical  and  are 

situated  in  a  horizontal  plane. 

b)  A  single  signal  is  present  and  has  the  form 

of  a  plane  wave  traveling  through  an  homogen¬ 
eous,  linear  and  nondispersive  medium. 
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c)  The  seismometer  outputs  are  identical  except  for 
a  time  delay  due  to  the  finite  phase  velocity 
of  the  wave. 

The  signal  useds  P  piase  of  L0NGSH0T  recorded  at  the 

center  element  of  the  LASa  A^  cluster. 

o 

The  patterns  were  corrected  to  include  the  effects  of 
the  attenuation  of  the  signal  during  the  path  from  the  hypo- 
center  to  the  array,  and  the  effect  of  the  geometric  spreading. 

Suppose  that  the  array  is  formed  by  N  seismometers  and 
let  their  locations  defined  by  the  vectors  r^,  k  =  1...,  N 

The  output  of  the  ktt  seismometer  will  be  % 

Xk(t)  =  S  (t  -  _a.'rk) 

where  S(t)  is  the  transient  used 

*a(t)  is  the  delay  vector  per  unit  distance 
Por  a  given  array  location,  every  point  of  the  earth  not 
situated  in  the  core-shadow  region  is  characterized  by  a  vector 
cT.  Let  7x  be  the  vector  associated  with  a  particular  point 
of  the  earth  that  we  shall  call  "Pocus"  of  the  array.  A  beam 
output  is  then  formed  by  delaying  all  seismometer  outputs  by 
a  •  Tt ,  k  =  1...,  N  and  summing  these  delayed  signals. 

The  Pocus  is  characterized  by  two  coordinates  with  respect 
to  the  center  of  the  arrays 

A0  =  Range  to  the  center  of  the  array 

Azq  =  Az  of  the  great  circle  that  joins  the 
Pocus  with  the  center  of  the  array 
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A  P-wave  signal  that  reaches  the  array  coming  from  Aq 
and  Azq  (The  Pocus )  will  produce  the  maximum  output  of  the 
array,  with  this  processing  scheme.  The  same  delays  (a  •  “r.  ) 
acting  on  a  signal  created  hy  an  event  of  the  same  magnitude 
hut  whose  hypocenter  is  not  located  at  the  focus  of  the  array 
will  produce  an  output  that  is  distorted  and  smaller  in 
amplitude# 

The  Radial  beam  pattern  computed  shows  how  the  maximum 
amplitude  of  the  delayed -sum  varies  as  the  event's  hypo- 

center  moves  away  from  the  focus,  over  the  great  circle 

As  expected,  effect  of  the  attenuation  and  the  geometric  spread¬ 
ing  acts  to  increase  the  sidelobes  at  ranges  smaller  than  the 

focus  and  their  effect  is  the  opposite  for  ranges  greater  than 

V 

The  transverse  pattern  computed  is  a  plot  of  the  maximum 
amplitude  of  the  delayed-sum  as  we  place  events  at  different 
points  on 
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One  output  of  the  program  is  the  transient  formed  by  the 
delayed  sum  of  the  seismometer  recordings  for  every  point  where 
the  beam  pattern  was  computed.  Figures  5 -a  to  5-e  show  these 
transients  for  ranges  varying  from  20°  to  96°  and  a  focal  range 
of  62°. 


Correction  for  Attenuation  and  Geometric  Spreading 

Geometric  Spreading: 


Ci  =  Correction  factor  =  — 
(amplitude)  ^ 


Aq  =  focus  range 


A  =  range 

This  correction  is  valid  for  a  spherically  symmetric 
traveling  wave. 

Attenuation  loss  correction; 

The  amplitude  of  seismic  signals  is  assumed  to  have  the 
following  range  and  frequency  dependence, Ref. 4  (after  removal  of 
geometric  spreading) ; 

U) 

A 

A  (uj,a)  =  e~ 

where : 


A  is  the  amplitude  of  the  wave. 

u)  is  the  angular  frequency  of  the  wave. 

A  is  the  distance  (central  angle)  between  source 
and  receiver. 

c  is  the  velocity  of  propagation. 

Q  is  the  ratio  of  energy  stored  to  energy  lost 
per  cycle  in  the  earth. 


Beam  Patterns  for  LAS A 

The  program  discussed  above  was  then  used  to  compute 
several  patterns  for  LASA  -  at  different  focal  ranges  and 
different  bearings  from  the  array.  For  these  computations,  it 
is  first  noted  that  LASA  is  formed  of  clusters  of  seismometers. 
The  dimensions  of  these  clusters  are  approximately  one  wave¬ 
length  at  the  predominant  frequency  of  teleseismic  signals . 
Thus,  the  pattern  for  this  array  can  be  found  by  multiplying 
the  pattern  obtained  using  each  cluster  as  a  single  element  of 
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an  array  with  the  patterns  of  a  single  cluster.  Since  the 
size  of  the  clusters  is  small  compared  to  the  wavelength  of  the 
signals  received,  the  beamwidth  of  these  subarrays  is  such  as 
to  include  all  teleseismic  arrival  angles.  Thus  the  pattern  for 
the  whole  array  can  differ  at  most  by  a  few  dB  from  that  computed 
for  the  center  elements  of  the  clusters  alone,  and  thus  the 
patterns  given  are  for  a  21  element  array  with  the  elements 
located  at  the  positions  of  the  center  of  the  clusters.  The 
following  three  foci  have  been  chosen  as  representative  ofLASA 
patterns  s 


1.  A0  =  62° 
Azq  =  277° 

2.  A0  =  56° 

Az  =  264° 
o 

3.  A0  =  64° 

Az  =  215° 
o 


Figures  6  and  7 


Figures  6  and  9 


Figures  10  and  11 


Plots  of  the  radial  beam  patterns  (Figures  6,6,  and  10) 
show  the  effect  of  the  attenuation  and  geometric  spreading  on 
the  pattern.  The  mainlobe  is  not  strongly  affected  but  the  side- 
lobes  change  appreciably. 

Plots  of  the  transverse  beam  patterns  (Figures  7  ,  9  , 11  ) 
of  LASA  center  elements  do  not  show  such  strong  effects. 

The  3-dB  beamwidths  obtained  are  the  following: 


with  attenuation  and 
geometric  spreading 


without  attenuation  and 
geometric  spreading 


A0  =  62° 
Azq  =  277° 

A  rr*0 

A0  =  5o 
AzQ  =  264° 


©r  =  13° 


0r  =  10.5° 

°t  =  5.7° 


©r  =10.0° 
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with  attenuation  and  without  attenuation  and 
geometric  spreading  geometric  spreading 

A0  =  64°  0r  =  13.75°  ©r  =  17,00° 

kzn  =  215°  ©+  =  12.5° 

Because  of  the  poor  sidelobes  of  the  LASA  patterns,  parti¬ 
cularly  in  the  radial  direction,  an  experiment  was  performed  to 
determine  the  degree  to  which  sidelobes  ould  be  controlled  if 
one  were  not  restricted  to  clustering  of  elements.  Such  an 
array  is  not  proposed  at  this  time  for  a  number  of  reasons,  how¬ 
ever,  the  patterns  it  could  produce  are  instructive. 

The  array  chosen  for  this  experiment  was  one  whose  elements 
(525)  were  distributed  randomly  (normal  distribution  in  both  x 
and  y  coordinates)  over  a  circle  of  radius  1°. 

Random  Distribution  of  Seismometers 

It  will  be  assumed  that  the  coordinates  which  locate  the 
array  elements  are  independent,  zero  mean,  normal  random  vari¬ 
ables  o  In  the  following,  the  density  function  of  a  zero  mean, 

2 

normal,  random  variable  with  variance  a  will  be  denoted  by 
N(o,a).  In  practice,  due  to  the  finite  size  of  the  array,  the 
elements  will  follow  a  truncated  normal  distribution  rather 
than  a  true  normal  distribution.  Two  parameters  define  a  sym¬ 
metric  truncated  zero  mean  normal  distribution: 

a2  =  variance 

=  maximum  value (trunca ti on  point) 

In  our  case  ^  corresponds  to  the  maximum  distance  from  the 
array  center  at  which  an  array  element  can  be  located. 

In  determining  the  array  pattern  for  a  radial  cut  at  an 
angle  ©  (see  Figure  below)  it  is  first  necessary  to  project  the 
array  elements  onto  the  radial  line.  Hence  the  projection  of 
the  array  element  with  coordinates  (x,y)  is  given  by 

R  =  xsin  ©  +  ycos  0 

If  x  and  y  are  independent,  random  normal,  variables  with 
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zero  means,  then  R  will  also  be  a  zero  mean,  normal,  random 
variable  since  it  is  a  linear  combination  of  normal  random 
variables.  Also,  becatise  of  the  sine-cosine  coeffecients,  the 
variance  of  R  will  be  the  same  as  those  of  x  and  y. 


Y 


Hence  for  normally  distributed  elements,  we  expect  array  patterns 
that  are  independent  of  Azq»  the  azimuth  of  the  focus. 

The  use  of  randomly  distributed  elements  as  discussed  above 
will  result  in  an  amplitude  taper  for  any  radial  cut  and  hence 
a  lower  sidelobe  level  (compared  to  the  clustered  array)  should 
be  obtained.  Figure  12  shows  such  a  distribution,  with  a  =0.5 
degrees  and  R^x  =  1  degree.  The  relation  between  a  and 
was  chosen  to  be  o  =  £  R^  in  order  to  avoid  strong  sidelobe 
generation  due  to  truncation  of  the  density  function.  For 
o  >  £  Rfiax  the  actual  array  beamwidth  will  be  smaller  than  com¬ 
puted,  but  the  truncation  effect  will  start  to  become  important. 
No  attempt  has  been  made  to  find  an  optimization  criterion. 
Figures  13  and  14  show  the  computed  transient  patterns  for  the 
random  array  of  Figure  11  when  focussed  at  a  range  of  56°  at  an 
azimuth  of  264°.  Also  shown  in  these  figures  are  the  comparable 
patterns  for  the  present  LASA  configuration. 
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Radial  Beam  Pattern  for  Randomly  Distributed  Elements 

In  analyzing  the  beam  pattern  of  a  random  LASA  array,  a 
number  of  simplifying  assumptions  will  be  made.  The  array  will 
be  treated  as  a  continuous  distribution  rather  than  discrete, 
truncation  effects  will  be  ignored,  and  the  incoming  wave  will 
be  assumed  to  be  monochromatic  rather  than  broadband*  With  the 
above  simplifications  it  will  be  shown  that  the  resulting  array 
pattern  of  the  random  array  is  essentially  the  Eourier  trans¬ 
form  of  the  probability  density  function  of  the  random  variable 
R.  Hence,  for  the  assumed  normal  distribution  of  array  elements, 
the  resulting  array  pattern  will  have  a  Gaussian  shape. 

As  discussed  previously  if  the  coordinates  x,  y  of  an  array 
element  are  independent  random  variables  with  normal  distributions 
N(o,o),  the  projection  R  of  the  element  R  on  a  radial  line  will 
also  have  a  normal  distribution  N(o,a)»  For  a  monochromatic 
signal  the  radial  pattern  cp(a)  is  given  by 


cpt(a) 

where 


1 

N 


N 

2 


i=i 


j2rrfaR1 


f  =  frequency  of  the  incoming  wave 
N  =  number  of  array  elements 

a  =  relative  delay  per  unit  distance  along  the  radial 
R^  =  the  projection  of  the  ith  element  onto  the  radial. 

It  should  be  noted  that  a  the  relative  delay  will  be  a  function 
of  9,  the  angle  of  the  radial,  and  the  phase  velocity  of  the 
incoming  signal*  Now  since  the  R^  are  random  variables,  it  is 
the  average  array  pattern  cp'a)  that  is  of  interest*  Hence 


cp(a)  = 


N  jZnfaR. 
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Thus : 


cp(a)  = 


|  J  ef2"faR  1  dR  1  = 

32nafR 

e 

fStTo 

Hence  cp(a)  is  given  hy  the  magnitude  of  the  characteristic 
function  of  the  probability  density  function  of  the  random 
variable  R.  The  characteristic  function  of  a  Gaussian  random 
variable  has  a  Gaussian  form.  Hence  2 


cp(a)  =  e 


0_2  2.2  2 
2tt  a  f  a 


where 


k  = 


2rr 


a  * 


Hence  cp(a)  is  equal  to  k  N  (o,  2rrof  ).  Due  to  the  reciprocal 
properties  of  the  Fourier  transform  the  beamwidth  of  the  random 
array  is  inversely  proportional  to  the  standard  deviation  of 
the  element  locations. 

For  the  expected  LASA  signals,  the  received  energy  will 
be  concentrated  in  the  vicinity  of  one  cycle.  Hence  for 
illustration,  the  array  pattern  has  been  calculated  on  the  basis 
of  a  received  signal  frequency  of  1  Hz  .  Results  of  the  compu¬ 
tation  (See  Figure  13)  show  that  the  assumptions  made  above  were 
valid. 
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ARRAY:  CENTER  ELEMENTS  OF  LASA  CLUSTERS 
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FIGURE  7  -  TRAVERSE  PATTERN 
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FIGURE  9  -  TRANSVERSE  PATTERN 
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FIGURE  11  -  TRANSVERSE  PATTERN 
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Section  II-C  DIMUS  Array  Processing 

Many  investigators  have  treated  the  case  of  clipped  (use 
of  zero  crossing  information  only)  processing  of  signals  in  the 
presence  of  noise.  Anderson  of  Scripps  Oceanographic  Institute 
coined  the  name  of  DIMUS  (digital,  multiple  steered)  array  pro¬ 
cessing  for  additive  arrays  using  clipper  amplifiers  and  shift 
registers  for  delay  compensation.  Anderson  concluded  that  the 
loss  in  array  processing  for  the  DIMUS  is  1.95  dB  while 

authors  have  shown  that  this  loss  may  he  higher  or  lower  depend¬ 
ing  on  the  character  of  the  interfering  noise  ^  .  These 

losses  in  processing  gain  are  for  the  case  of  low  signal-to- 
noise  ratio  at  each  of  the  individual  sensors  which  make  up  the 
array.  Eor  higher  input  signal-to-noise  ratios,  this  loss  can 
become  much  more  severe.  This  may  readily  be  seen  from  the  fact 
that  the  DIMUS  output  is  constrained  to  lie  somewhere  between 
positive  and  negative  limits  equal  to  the  number  of  sensors  in 
the  array.  Where  the  problem  is  the  detection  of  weak  signals 
in  noise,  a  loss  in  output  signal-to-noise  ratio  of  some  two  dB 
or  less  in  considered  insignificant  for  seismic  work.  This  is 
because  of  the  tremendous  variability  of  seismic  measurements. 
What  does  appear  significant,  however,  is  the  distortion  intro¬ 
duced  for  large  input  signal-to-noise  ratios  (S/N>1) .  In  this 
case,  discrimination  between  earthquakes  and  underground  ex¬ 
plosions  on  the  basis  of  the  received  waveforms  might  be  severly 
curtailed.  Hence,  a  study  of  DIMUS  processing  must  be  under¬ 
taken  from  two  distinct  points  of  view;  detection,  and  identi¬ 
fication. 

Before  proceeding  further  with  a  discussion  of  DIMUS  pro¬ 
cessing  itself,  let  us  first  mention  several  reasons  why  such 
processing  might  have  distinct  advantages  for  a  large  aperture 
seismic  array.  Most  obvious  is  the  reduction  in  the  amount  of 
data  which  must  be  handled  and/or  stored.  Since  in  DIMUS  each 
sample  conveys  one  bit  of  information,  there  is  available  an 
immediate  reduction  of  some  15  to  1  in  the  data  handled  over 
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that  necessary  in  the  current  design  of  the  Montana  1ASA  system. 
.Further,  if  we  might  limit  our  attention  to  seismic  frequencies 
below  two  cycles  per  second,  the  sampling  rate  itself  could  be 
reduced  to  some  four  or  five  samples  per  second  rather  than  the 
twenty  samples  per  second  for  LASA  Montana.  The  net  result  of 
such  changes  would  be  a  data  rate  some  sixty  times  less  than 
current  practice.  One  result  of  such  a  reduction  in  data  rate 
would  be  the  possibility  of  storing  the  entire  output  of  a  LASA 
on  only  two  ten  inch  recording  tapes  each  day.  Perhaps  more 
important,  for  a  500  seismometer  array,  the  total  output  data 
rate  (taking  each  seismometer  separately)  would  be  approximately 
2.5  kilobits  per  second.  Such  a  rate  is  clearly  low  enough  to 
permit  the  transmission  of  the  total  raw  data  from  such  a  site 
to  a  central  data  processing  facility  over  a  single  telephone 
circuit  in  real  time.  Similar  dramatic  reductions  in  beamform¬ 
ing  equipment  are  also  easy  to  propose.  Finally,  if  such  a 
drastic  reduction  in  data  rate  from  individual  instruments 
proves  completely  satisfactory,  it  may  be  realistic  to  con¬ 
sider  arrays  where  the  seismometers  need  not  be  clumped.  Rather, 
each  instrument  could  be  self  powered  by  a  fuel  cell  or  solar 
battery,  and  the  instrument  output  relayed  to  the  array  center 
by  means  of  very  simple  communications  equipments  located  at 
the  instrument  sites.  So  much  for  the  potential  advantages  of 
DIMUS  processing  itself.  As  stated  above,  we  will  consider 
such  processing  from  two  standpoints;  detection  and  identifica¬ 
tion. 

Detection  Aspects  of  DIMUS  Beamforming  for  LASA 
As  stated  earlier  in  this  section,  previous  studies  of  the 
detection  capabilities  of  DIMUS  beamforming  systems  have  pre¬ 
dicted  a  one  to  two  dB  loss  in  processing  gain  for  such  systems 
when  compared  to  analog  systems,  depending  on  the  characteristics 
of  the  interfering  noise.  While  microseismic  noise  is  admittedly 
non-white  and  non-stationary  in  the  long  run,  it  has  been  hypoth¬ 
esized  that  at  least  for  short  intervals  of  time,  this  noise  is 
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Gaussian  (see  Section  II-D  for  example).  For  this  kind  of  noise, 
one  would  expect  that  the  DIMUS  processing  would  be  two  dB  poorer 
than  analog  processing  given  low  input  signal-to-noise  ratios. 

To  test  this  hypothesis,  a  continental  size  seismic  array  was 
formed  using  the  five  continental  U.S.  array  stations  and  NPNT 
as  clusters.  Each  of  the  clusters  was  summed  without  delay, 
and  the  phased  sum  taken  of  these  clusters.  The  event  taken  for 
the  experiment  was  a  surface  focus  event  of  magnitude  4.6,  with 
epicenter  located  in  Algeria.  Differential  delays  for  the  clus¬ 
ters  were  determined  from  another,  larger,  surface  focus  event 
located  in  the  same  general  vicinity. 

Figure  15  is  the  record  from  a  single  seismometer  located 
in  the  WMO  "cluster".  This  record  is  typical  of  those  for  each 
instrument  in  the  large  array.  From  this  figure,  it  is  easy  to 
see  that  the  criterion  of  low  signal-to-noise  ratio  at  each 
seismometer  was  met  for  this  test.  Figure  16  shows  both  the 
analog  sum  of  the  seismometers  for  the  array  focussed  at  the 
event  under  study,  the  sum  of  the  hardlimited  seismometer  out¬ 
puts,  and  also  the  same  two  signals  after  band  pass  filtering 
(0.75  to  1.5  cps).  This  filter  served  the  dual  purposes  of 
removing  out  of  band  interference  and  harmonic  distortion  intro¬ 
duced  by  the  limiting  process.  The  event  in  question  is  clearly 
visable  on  each  of  these  four  traces,  however,  the  final  signal- 
to-noise  ratio  leaves  something  to  be  desired. 

Since  the  peak  outputs  of  both  kinds  of  array  processing 
have  been  normalized  in  the  traces  shown,  it  was  easy  to  com¬ 
pute  the  rms  value  of  the  noise  for  both  array  outputs  immedi¬ 
ately  preceding  the  arrival  of  tne  event.  The  ratio  of  clipped 
noise  to  unclipped  noise  was  found  to  be  2.1  or  about  the  ex¬ 
pected  loss  in  DIMUS  processing.  While  the  length  of  the  noise 
record  available  for  analysis  on  each  of  these  records  was  small 
(approximately  15  seconds),  and  only  one  experiment  was  per¬ 
formed,  the  results  agree  so  well  with  theoretical  expectations 
that  one  may  assume  that  the  detection  performance  of  DIMUS 
beamformers  can  be  predicted. 
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Before  leaving  the  detection  aspects  of  the  DIMUS,  it  is 
interesting  to  note  that  such  beamformers  make  possible  CFAR 
(constant  false  alarm  rate)  automatic  detectors  without  tht 
usual  problems  of  automatic  gain  controlled  amplifiers.  To 
see  how  this  is  so,  first  consider  that  the  noise  imputs  to 
eacn  of  the  seismometers  are  independent  and  have  zero  means. 
Then,  the  sum  of  N  hard  limited  instrument  outputs  will  range 
from  -N  to  +N  and  possesses  a  binominal  probability  distribution 
with  zerc  nean  -  no  matter  what  the  strength  of  the  noise  at 
each  of  the  seismometers.  It  is  easy  to  see,  then,  that  a 
constant  threshold  set  at  the  array  output  will  give  a  con¬ 
stant  and  predictable  false  alarm  rate. 

To  determine  the  probability  of  detection  for  such  a  system, 
first  consider  a  constant  signal  at  the  input  to  each  seismo¬ 
meter.  Then  in  the  absence  of  noise,  all  of  the  hard  limiters 
following  the  seismometers  would  have  the  same  polarity,  and 
the  magnitude  of  the  output  of  the  summer  of  the  system  would 
be  exactly  N.  When  sufficient  noise  is  present  at  the  inputs 
to  the  seismometers,  not  all  of  the  hard  limited  outputs  will 
be  of  the  same  polarity.  Here  the  magnitude  of  the  output  of 
the  array  summer  can  take  values  less  than  N,  the  exact  value 
being  a  function  of  the  signal-to-noise  ratio  at  the  inputs  to 
the  individual  seismometers. 

Figure  17  illustrates  how  it  is  possible  to  determine  the 
probability  distribution  of  the  magnitude  of  the  summer  output 
as  a  function  of  the  input  signal-to-noise  ratio.  The  figure 
shows  the  distribution  of  signal  plus  noise  at  the  input  to  a 
single  sensor  as  a  normal  distribution  with  mean  value  E&, 
the  signal  voltage.  This  results  from  the  assumption  that  the 
noise  is  normally  distributed  with  zero  mean.  From  the  figure, 
it  is  easy  to  see  that  the  probability  that  noise  will  cause 
the  limiter  output  to  be  of  different  polarity  than  would  be 
true  in  the  noiseless  case  is  simply  the  area  of  the  distribu¬ 
tion  for  negative  values  of  the  variable.  Thus,  given  the 
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input  signal-to-noise  ratio  at  the  individual  seismometers, 
it  is  possible  to  compute  the  probability  of  noise  producing  an 
incorrect  polarity  for  a  single  seismometer  output.  One  can  in 
fact  compute  what  this  signax-to-noise  ratio  must  be  for  a  given 
probability  that  no  more  than  m  of  the  N  outputs  will  be  so 
reversed.  If  the  threshold  is  set  at  a  value  N-m,  this  proba¬ 
bility  will  be  the  probability  of  detection.  Figure  16  give 
several  illustrations  of  the  detection  performance  (probability 
of  detection  versus  input  signal-to-noise  ratio  for  a  given 
false  alarm  rate)  for  a  DIMUS  CFAR  detector.  The  illustrations 
are  for  a  21  element  array  and  a  49  element  array,  the  first 
chosen  to  illustrate  what  1ASA  performance  would  be  if  hard 
limiting  takes  place  after  analog  summation  of  the  individual 
seismometers  of  the  clusters.  As  can  be  seen  from  the  figure, 

the  21  element  CFAR  possess  too  large  a  false  alarm  rate  for 

most  purposes.  Increasing  the  number  of  elements  by  a  factor 
of  approximately  two  improves  the  false  alarm  rate  by  two 
orders  of  magnitude  to  a  value  more  consistant  with  good  system 
operation  (Pfa  **  10”^). 

Identification  With  A  DIMUS  Array  Processor 
As  stated  at  the  beginning  of  this  section,  the  DIMUS 
beamformer  is  inherently  a  nonlinear  device.  A  recent  paper 
by  Remley^^,  presents  the  first  derivation  of  the  transfer 
function  for  such  a  processor.  This  transfer  function  is  given 

in  Figure  19,  where  the  output  signal  voltage  is  given  as  a 

function  of  the  input  signal-to-noise  ratio  at  each  limiter. 

As  can  be  seen  from  the  figure,  when  the  signal-to-noise  ratio 
is  on  the  order  of  unity  or  less,  the  output  signal  waveform 
is  expected  to  be  linearly  related  to  the  input  wave-form. 

Where  the  number  of  elements  in  the  array  is  large,  not  only  will 
the  output  signal  be  a  faithful  representation  of  the  input 
signal,  but  the  output  signal-to-nosie  ratio  may  also  be  large  - 
even  for  low  input  signal-to-noise  ratios.  Here,  the  DIMUS 
output  may  be  used  exactly  as  the  analog  output  is  in  the  de¬ 
termination  of  the  source  of  the  seismic  signals.  Of  course 
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the  signal-to-noise  ratio  for  the  DIMUS  will  he  somewhat  less 
(2  dB)  than  for  the  conventional  array,  however,  no  distortion 
of  the  waveshape  will  have  been  introduced  by  the  nonlinear 
procef  sing. 

When  the  input  signal-to-noise  ratio  i3  high,  however,  con¬ 
siderable  distortion  of  the  seismic  waveform  may  take  place. 

This  would  be  true  for  magnitude  5  events,  for  example.  To 
discover  the  magnitude  of  such  effects,  the  records  of  some 
13  teleseismic  events  (both  earthquakes  and  underground  explo¬ 
sions)  were  processed  with  both  conventional  and  DIMUS  techniques. 
At  least  10  widely  spaced  seismometers  were  used  to  form  beams 
focused  at  each  event.  The  signal-to-noise  ratio  was  high  in 
all  cases.  Thus  the  differential  time  dei^rs  could  be  picked  by 
eye  from  the  individual  seismograms.  The  measure  chosen  to 
determine  the  similarity  between  the  analog  and  DIMUS  outputs 
was  the  zero  offset  correlation  coefficient  for  the  P-coda 
signals  of  these  events.  Here  P-coda  was  taken  to  be  that 
portion  of  the  seismogram  from  three  to  thirteen  seconds  follow¬ 
ing  the  onset  of  the  P-phase.  The  choice  of  this  measure  was 
based  on  the  hypothesis  that  the  identification  of  source  mechan¬ 
ism  would  be  based  mainly  on  the  waveforms  found  in  this  time 
interval  of  the  seismogram,  although  it  was  recognized  that  most 
complexity  criteria  employed  to  date  required  a  measurement  of 
the  ratio  of  main  P  energy  to  energy  in  specific  portions  of  the 
P-coda.  Undoubtably  this  choice  of  hypothesis  was  strongly 
affected  by  the  interesting  results  of  study  of  the  P-coda 
correlation  coefficient  as  a  discriminant  (see  Section  III). 

Results  of  these  computations  were  much  as  expected.  P-coda 
was  emphasized  with  respect  to  main  P,  and  the  general  noise 
level  increased.  These  two  effects  are  illustrated  by  Figure  20, 
which  shows  the  filtered  (0.75  to  1.50  cps)  analog  and  DIMUS 
beam  outputs  (14  seismometers)  for  a  typical  sarthquake.  What 
is  somewhat  unexpected  is  the  remarkable  similarity  of  the  coda 
waveforms  for  the  two  signals.  In  fact,  the  coda  correlation 
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coefficient  is  0,96.  Of  the  13  events  studied,  this  coda  corre¬ 
lation  coefficient  ranged  from  0.97  to  0.59,  with  only  two  events 
showing  correlation  coefficients  lower  than  0.76.  With  this  high 
correlation,  it  is  hypothesized  that  the  coda  correlation  dis¬ 
criminant  discussed  in  Section  III  could  he  satisfactorily  em¬ 
ployed  with  convent? onal  DIMUS  processing.  There  are,  however, 
a  number  of  ways  by  which  the  nonlinear  transfer  function  of  the 
conventional  DIMUS  may  be  linearized. 

To  date,  one  test  of  the  simplest  modification  of  DIMUS  pro¬ 
cessing  has  been  made.  Recognizing  that  only  events  with  large 
signal-to-noise  ratio  will  be  distorted  by  the  conventi  onal 
DIMUS,  it  was  hypothesized  that  the  output  of  a  single  seismo¬ 
meter  of  the  array  should  be  a  good  estimate  of  the  desired 
output  signal  -  particularly  for  the  largest  portions  of  the 
record.  Thus,  in  addition  to  the  hard  limited  outputs  of  the 
individual  seismometers,  one  analog  channel  was  used  to  estimate 
the  output  of  a  conventional  analog  array.  The  modified  array 
output  was  then  computed  in  the  following  manner.  At  the  first 
point  in  time  that  all  hard  limited  signals  potsess  the  same 
polarity  (presumably  that  time  for  which  the  input  signal  entered 
the  region  of  severe  distortion  in  the  transfer  function),  the 
corresponding  value  of  the  single  analog  channel  was  noted.  The 
analog  channel  was  then  normalised  so  that  its  value  at  the  noted 
point  was  equal  to  the  number  of  seismometers  making  up  the  array 
(21  for  the  case  studied).  Then,  for  those  values  of  the  sum 
of  the  haid  limited  channels  was  equal  to  the  number  of  channels, 
the  hard  limited  sum  was  replaced  by  the  corresponding  value  of 
the  normalized  analog  channel. 

Results  of  this  simple  "fix"  are  illustrated  in  Figure  21. 
Shown  in  the  figure  are  the  analog  beam  output,  the  modified 
DIMUS  beam  output,  and  the  conventional  DIMUS  bean  output.  From 
the  figure,  it  is  easy  to  see  that  even  this  simple  "fix" 
eliminated  nearly  all  of  the  distortion  introduced  by  straight 
DIMUS  processing.  Some  minor  distortion  still  exits  for  the 
main  P  arrival  and  the  first  second  of  the  coda,  however,  the 
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remaining  portion  of  the  modified  array  output  is  as  nearly 
identical  with  the  corresponding  portion  of  the  analog  processor 
as  could  be  desired.  In  addition,  it  can  be  seen  from  the  fig¬ 
ure  that  the  level  of  the  noise  in  the  modified  record  is  re¬ 
duced  over  that  for  the  conventional  DIMUS.  This  is  to  be 
expected,  as  it  will  be  the  ratio  of  low  energy  coda  to  noise 
which  is  expected  to  remain  unchanged  by  this  processing  tech¬ 
nique.  As  a  result,  when  the  coda  level  is  reduced  to  its 
proper  relationship  with  respect  to  main  P,  the  noise  should 
also  shrink  witti  respect  to  P.  Other,  more  sophistocated  modi¬ 
fications  of  the  DIMUS  technique  are  currently  under  study  with 
the  view  of  correction  of  the  distortion  in  the  first  three 
seconds  noted  above.  Even  without  such  further  modifications, 
the  work  to  date  strongly  suggests  that  the  inclusion  of  a 
single  analog  channel  with  a  DIMUS  processor  can  provide  all  of 
the  information  possible  with  an  analog  array,  and  it  can  do  so 
with  far  greater  simplicity. 
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FIGURE  15  -  ALGERIAN  EVENT  -  SINGLE  SEISMOMETER  WMO 


ANALOG  UNFILTERED 
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FIGURE  16  -  ANALOG  AND  DIMUS  BEAMS 
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FIGURE  17  -  PROBABILITY  DISTRIBUTION  FOR  SIGNAL  PLUS  NOISE 
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FIGURE  18  -  REPRESENTATIVE  DIMUS  DETECTION  CAPABILITY 
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FIGURE  19  -  STATISTICAL  EQUIVALENT  OF  DIMUS  FOR  GAUSSIAN  NOISE 

WITH  RMS  o 


ANALOG  ARRAY 
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FIGURE  20  -  COMPARISON  OF  ANALOG  AND  DIMUS  PROCESSING 


ANALOG  BEAM 
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FIGURE  21  -  COMPARISON  OF  BEAM  PROCESSING  TECHNIQUES 
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SECTION  II-D.  Spatial  Correlations  Observed  Across  LASA  from  LONGSHOT 

Section  III  presents  preliminary  measure 
ments  on  P-coda  correlation  between  seismograms  spaced  1000  to 
7000  km  for  a  small  group  of  earthquakes  (12)  and  nuclear  explo¬ 
sions  (11).  Although  the  number  of  events  was  much  too  small  to 
form  conclusions,  the  data  separated  in  such  a  way  as  to  suggest 
that  P-coda  correlation  measurements  may  possibly  form  the  basis 
for  a  powerful  classification  procedure.  Many  more  events  are 
now  being  prepared  to  test  the  hypothesis  of  that  note. 

(Briefly  stated,  the  mean  value  of  the  cross-correlation 
of  P-coda  for  the  explosions  fell  close  to  zero,  grouped  around 
0.3  to  0.4  for  shallow  earthquakes( defined  here  as  above  the  Moho) 
and  fell  in  between  for  earthquakes  deeper  than  the  Moho.  Since 
it  is  the  shallow  earthquake  which  it  is  most  difficult  to  separate 
from  a  surface  event,  the  large  separation  between  the  measured 
parameter  for  the  two  classes,  i.e.,  explosions  and  shallow  earth¬ 
quakes,  is  very  encouraging.) 

It  has  been  suggested  by  Mr.  L.  H.  Orpin  of  MITRE  that 
a  similar  experiment  be  conducted  over  a  region  of  LASA  dimensions 
(200  km),  as  the  results  may  conceivably  bear  on  LASA  array  design. 

This  work  has  been  started  and  the  first  numerical  results  are 
presented  here.  The  LASA  array  in  Montana  has  been  used 

as  the  test  bed  and  LONGSHOT  used  as  the  source.  The  distance 
from  the  center  of  the  array  to  the  source  was  about  5250  km. 

In  addition  to  P-coda  correlation  measurements  the 
spatial  cross-correlation  of  the  main  P  wave  and  of  noise  just 
prior  to  the  event  also  have  been  measured.  The  latter  quantities 
are  pertinent  to  prediction  of  array  processing  gain  relative  to 
the  detection  of  the  P  wave  in  noise.  The  coda  portion  used  in 
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the  experiment  was  3  to  13  seconds  following  P-onset,  typically 
after  main  P  died  out  but  before  pP  from  a  source  below  the 
Moho  would  arrive.  The  section  of  P  wave  used  in  the  calculation 
was  1  second  in  duration  centered  at  the  first  well-defined 
zero  crossing  of  the  P  wave.  Noise  was  10  seconds  in  duration 
from  12  seconds  to  2  seconds  prior  to  main  P. 

Our  cross-correlation  program  can  accommodate  20  seis¬ 
mograms.  We  selected  the  central  seismogram  of  each  cluster 
except  the  central  cluster.  Thus  each  of  the  four  subarrays  in 
the  B,  C,  D,  E,  and  F  rings  contributed  one  sensor.  The  number 
of  different  cross-correlation  coefficients  computable  from  N 
instruments  is  (N^-N)/2,  which  is  190  for  N=20.  The  results 
follow . 

P-Coda  Correlation 

The  mean  and  standard  deviation  of  the  190  correlation 
measurements  are  0.24  and  0.25.  The  distribution  of  P-coda  corre¬ 
lation  clearly  exhibits  distance  dependence,  however.  Figure  22 
shows  the  mean  values  in  the  10  km  intervals.  Table  I  lists 
the  means,  standard  deviations  and  the  number  of  correlation 
coefficients  in  each  interval. 

Although  the  variation  with  separation  between  the 
elements  is  erratic,  a  general  diminution  with  separation  is 
apparent.  The  mean  correlation  drops  toward  zero  in  the  neigh¬ 
borhood  of  100-150  km.  The  separation  or  distance  dependence 
becomes  clearer  in  Figure  23,  in  which  the  range  scale  is  made 
logarithmic.  Here  the  same  data  are  replotted  with  the  distance 
interval  one  decade  in  size.  The  decrease  in  the  erratic  nature 
(compare  with  Figure  22)  is  due  to  the  larger  number  of  sample 
points  in  each  interval. 
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TABLE  I 


MEANS,  STANDARD  DEVIATIONS 

,  AND 

NUMBER  OF  SAMPLES 

FOR  FIGURES  1 

AND  2 

DISTANCE 

INTERVAL 

P. 

a 

N 

(km) 

c 

8- 

.^7 

.26 

2 

10- 

•  37 

.24 

18 

20- 

•  39 

.21 

24 

30- 

.40 

.15 

20 

40- 

.29 

.24 

14 

50- 

.  26 

.16 

14 

60- 

.25 

.19 

16 

70- 

.24 

.20 

16 

80- 

.04 

.19 

11 

90- 

.05 

•  33 

14 

100- 

.13 

.20 

9 

110- 

.12 

.19 

10 

120- 

.01 

.04 

3 

130- 

.03 

.21 

4 

140- 

.16 

.17 

5 

150- 

-.13 

.27 

3 

l60- 

.10 

— 

1 

8-15 

.40 

.26 

11 

16-31 

•  37 

.19 

35 

32-63 

•  31 

.19 

52 

64-127 

.14 

.24 

76 

28-210 

.03 

.22 

14 
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On  the  assumptions  that : 

i.  the  coda  correlation  test  suggested  in  Section  III 
may  prove  useful,  and 

ii.  the  distance  dependence  measured  on  LONGSHOT 
is  somewhat  typical, 

it  is  pertinent  to  examine  the  distribution  of  subarrays  in  LASA 
so  as  to  estimate  the  number  of  instruments  which  are  adequately 
spaced  to  be  useful  for  the  test.  Such  examination  shows  that  the 
Aq  subarray  and  the  four  subarrays  in  the  P-ring  are  mutually  far 
enough  part.  Thus  from  a  single  LASA  a  five-station  test  could  be 
expected.  Our  experience  last  year  in  measuring  the  mean  coda 
correlation  indicates  that  about  twice  this  number  is  wanted.  Hence 
it  is  not  likely  that  a  single  LASA  would  suffice;  two  would  be 
needed  to  conduct  the  test  with  reasonable  confidence,  or  one  LASA 
plus  five  other  array  stations,  such  as  the  U.S.  or  Canadian  ob¬ 
servatories  . 

It  is  planned  to  test  the  latter  possibility. 

Signal 

The  cross-correlation  of  main-P  is  large.  The  mean  value 

is  0.97  and  the  standard  deviation  is  0.04.  The  large  correlation 

is  not  surprising  in  view  of  the  fact  that  main  P  correlation  was 

very  large  in  the  earlier  reported  work  where  the  distances  between 

9  10 

elements  were  thousands  of  kilometers.  * 

No  distance  dependence  is  observed  in  the  data  contrary  to 
the  expected  tendency  for  correlation  to  fall  off  with  separation 
between  seismograms.  The  smallest  cross-correlation  coefficient, 
which  is  0.72,  is  obtained  from  seismometers  in  C^  and  L^,  which 
are  only  about  40  to  45  kilometers  part.  All  C^  values  are  low 
compared  to  the  rest  of  the  cross-correlation  matrix.  The  mean 

correlation  of  C,  with  all  the  rest  is  0.87,  possibly  due  to  error 
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in  picking  event  arrival  time.  Dropping  on  this  assumption 

raises  the  mean  to  0.9S  and  reduces  the  standard  deviation  to  0.02. 

The  distribution  of  the  correlation  coefficient  (minus  C_) 

5 

is  plotted  in  Figure  24. 

Noise 

In  the  earlier  work  where  the  distances  between  seismo¬ 
meters  were  thousands  of  kilometers,  the  mean  noise  cross* 
correlation  appeared  unbiased  and  reasonably  symmetrical.  It 
was  of  interest  to  determine  if  the  same  was  true  when  the  distance 
between  elements  shrank  drastically.  Across  LASA  the  minimum 
spacing  is  about  S  km.  The  noise  on  the  seismic  records  just 
prior  to  the  10NGSH0T  event  permitted  such  a  test.  Figure  25 
shows  the  cumulative  distribution  plotted  on  probability  paper 
(dashed  curve).  The  mean  and  standard  deviation  are  0.009  and 
0.35. 

A  test  for  independence  of  the  noise  between  elements 
may  be  made  by  comparing  this  distribution  with  the  theoretical 
cumulative  distribution  for  the  cross-correlation  coefficient 
between  two  independent  Gaussian  sources  when  the  sample  size 
is  comparable  to  the  estimated  number  of  independent  samples  in 
the  noise  records  from  the  seismograms.  This  quantity  is 

N  =  2TW 

where  T  is  the  length  of  the  noise  record  and  ¥  is  the  noise 

bandwidth.  The  length  of  the  record  T  =  10  seconds.  From 

11 

other  spectral  measurements  of  other  records j 
1/4  <  ¥  <  1  cps. 

Thus  the  estimated  number  of  independent  noise  samples  in  each 
record  is  about  5  to  20  samples. 
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Shown  in  Figure  25  are  the  theoretical  curves  for  the 

cross-correlation  coefficient  between  two  independent  Gaussian 

sources  for  a  range  of  sample  size  which  brackets  the  estimated 

12 

number  of  independent  samples  in  each  record.  It  is  apparent 
that  the  measured  distribution  is  consistent  with  N  =  9,  except 
at  the  tail,  where  it  deviates  markedly.  It  is  believed  that 
this  wild  departure  is  attributable  to  the  limited  number  of 
correlation  coefficients  from  which  the  distribution  was  made 
up.  That  is,  the  experimental  distribution  had  a  very  small 
number  of  extreme  values.  If  this  is  true,  a  more  extensive 
test  will  eliminate  the  seemingly  abrupt  truncation. 

These  results  are  consistent  with  the  hypothesis  that 
the  noise  at  each  cluster  is  independent  of  that  of  every  other 
cluster . 

Comment  on  Detection  Efficiency 

It  was  indicated  in  Reference  2  that  the  gain  in  signal- 

to-noise  power  ratio  for  a  teleseismic  array  is  a  first  order 

function  of  three  quantities.  These  are  the  number  of  elements  N, 

the  mean  spatial  correlation  of  signal  p  ,  and  the  mean  spatial 

cross-correlation  of  noise  p  .  The  relation  is 

_  n 

p  (N-l)+l 

G  =  ^ - 

Pn(N-l)+l 

In  this  expression  N  is  the  number  of  clusters  or  subarrays  rather 
than  the  total  number  of  seismometers.  This  is  because  the  di¬ 
mension  of  each  subarray  is  only  the  order  of  a  wavelength. 

The  ideal  asymptotic  enhancement  is  obtainable  when 

p"  =1  and  p"  =0.  The  limiting  gain  in  this  case  is  N,  the 
s  n 

number  of  subarray  clusters. 


GENERAL  ATRONICS  CORPORATION 


It  is  to  be  expected  that  as  the  separation  between 
elements  increases,  both  mean  correlations  p  and  p  will  decrease. 
Thus  to  maximize  detection  efficiency  an  optimum  array  size  is  one 
which  is  not  so  large  that  p  departs  greatly  from  unity,  while  it 
is  large  enough  so  that  p"  is  near  zero.  The  LONGSHOT  experiment 
reported  above  suggests  that  the  LASA  dimension  may  be  large  enough 
to  accomplish  the  latter  while  not  so  large  as  to  cause  significant 
degradation  in  the  signal  correlation.  Thus  for  detection  of  main 
P  in  microseismic  noise,  200  kilometers  may  well  be  within  the  broad 
optimum  array  size  range. 
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FIGURE. ,22  -  SEPARATION  BETWEEN  ELEMENTS  (KM) 


FIGURE ...23  -  SEPARATION  BETWEEN  ELEMENTS  (KM) 

MEAN  SPATIAL  CORRELATION  OF  LONGSHOT  CODA  AS  A  FUNCTION 
OF  THE  DISTANCE  BETWEEN  SEISMOMETERS.  LENGTH  OF  CODA 
IS  10  SECONDS  BEGINNING  3  SECOND  AFTER  P-ONSET 
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FIGURE  .]2H  -  DISTRIBUTION  OF  THE  SPATIAL  CORRELATION  ACROSS 
LASA  OF  FIRST  SECOND  OF  LONGSHOT  P-WAVE 
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FIGURE  25  -  CUMULATIVE  DISTRIBUTION  OF  THE  SPATIAL  CROSS¬ 
CORRELATION  ACROSS  LASA  OF  10  SECONDS  OF  NOISE 
PRIOR  TO  LONGSHOT.  SOLID  CURVES  ARE  THEORETICAL 
DISTRIBUTIONS  FOR  INDEPENDENT  NORMAL  SOURCES. 

N  =  NUMBER  OF  INDEPENDENT  SAMPLES  IN  EACH  RECORD 
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Section  III  Coda  Correlation  as  a  Discriminant 

Eor  several  months  General  Atronics  has  been  studying  the 
average  cross-correlation  coefficient  between  seismograms 
from  widely  separated  stations  in  response  to  common  events. 

This  has  been  done  for  microseismic  noise,  P  wave  onset, 

P-coda,  and,  where  possible,  pP^  All  seismograms  were  time 
aligned  by  eye.  The  separation  between  stations  ran  up  to 
about  7000  kilometers. 

The  purpose  of  forming  these  mean  correlation  coefficients 
was  to  determine  the  extent  to  which  signal  enhancement  would 
be  achieved  by  team  formation  from  arrays  of  enormous  extent. 

Por  the  detection  problem  the  signal  was  considered  to  be  the 
P  wave  in  the  neighborhood  of  onset  and  the  noise  was  considered 
to  be  microseismic?  for  the  classification  problem  the  signal 
was  considered  to  be  pP  and  the  noise  or  background  interference 
was  considered  to  be  P-coda. 

Array  theory  indicates  that  it  is  the  mean  cross-correla¬ 
tion  between  appropriate  portions  of  the  seismogram  which  domin¬ 
ates  the  amount  of  signal  enhancement  relative  to  background 
interfeience  provided  by  a  large  array.  Consequently,  the  mean 
values  of  the  four  cross-correlation  coefficients  were  obtained. 

Still  speaking  historically,  it  was  found  that  the  mean 
cross-correlation  coefficients  for  the  P  wave  were  reasonably 
close  to  unity,  while  those  for  microseismic  noise  were  very 
close  to  zero.  Consequently,  excellent  signal-to-noise  enhance¬ 
ments  could  be  obtained  for  the  detection  problem.  Por  the 
classification  problem,  however,  the  data  were  much  less  con¬ 
clusive.  It  appeared  that  the  correlation  for  P-coda  between 
stations  separated  by  thousands  of  kilometers  was  much  too  high 
to  achieve  useful  enhancement  of  pP  in  P-coda.  However,  the 
data  were  too  fragmentary  to  form  such  a  firm  conclusion.  Con¬ 
sequently,  more  events  are  now  being  prepared  for  inclusion  in 
the  study. 
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Even  though  the  number  of  events  was  not  large,  the 
larger-tnan-expected  ^°da  correlation  demanded  examination 
of  the  data  in  any  wcys  which  night  disclose  consistencies, 
peculiarities,  etc*  One  such  examination  did  disclose  a 
potentially  important  discriminator  between  shallow  earthquakes 
and  surface  events „  This  is  the  subject  of  this  section  of  the 
report. 

A  total  of  25  events,  12  earthquakes  and  11  surface  events 
have  been  studied  to  date.  The  stations  were  widely  distributed 
and  were  in  no  sense  special.  The  section  of  coda  used  is  3  to 
13  seconds  after  P  onset.  This  interval  is  beyond  the  main  P 
wave  but  before  interaction  with  the  Moho  occurs.  It  thus  in¬ 
cludes  pP  for  the  range  of  earthquake  deptns  in  which  we  are 
most  interested. 

The  10  second  P-coda  interval  contains  200  samples  per 
record.  An  average  of  about  10  stations  each  recorded  each 
event;  hence  each  event  is  represented  by  about  2000  data 
samples.  Each  set  of  10  records  provided  45  f(N  -N)/2  correla¬ 
tion  coefficients.  The  average  of  these  45  numbers  is  the 
statistical  variate  of  interest.  It  is  designated  “  • 

The  total  number  of  samples  is  deceptively  large,  being 
about  2000  samples  for  each  of  20  events,  or  40,000  samples. 

More  important  to  remember  is  that  they  represent  only  20  events, 
which  is  an  inadequate  number  of  permit  conclusions  to  be  drawn. 
The  data,  however,  force  the  creation  of  a  hypothesis  which  is 
worthy  of  testing,  and  for  which  testing  is  now  underway. 

Figure  26  is  a  plot  of  the  fraction  of  events  having  a 
mean  coda  correlation  less  than  the  abscissa.  It  is  plotted 
for  each  of  the  two  classes.  The  points  for  the  23  events  are 
shown.  Also  shown  on  the  figure  are  the  general  locations  of 
the  surface  events,  the  depth  of  the  earthquakes,  and  the  mag¬ 
nitudes  for  both  classes  of  events.  As  can  be  seen  from  the 
figure,  no  particular  location  dependence  can  be  observed  in  the 
data.  Neither  can  a  magnitude  dependence  be  discovered.  There 
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can  be  observed,  however,  a  definite  relationship  between 
the  earthquakes  of  differing  depth.  More  will  be  said  about 
this  later. 

The  coordinates  for  Eigure  26  are  those  of  probability 
paper.  A  Gaussian  distribution  plots  as  a  straight  line  on 
this  paper.  The  number  of  events  is  very  small  in  each  of  the 
two  cases;  however,  if  we  assume  that  these  small  numbers  of 
event s  typify  the  real  world  distributions,  a  powerful  discrim¬ 

inant  test  becomes  available.  This  is  better  seen  in  Eigures  27 
and  28o  The  curves  of  Figure  27  are  the  derivatives  of  the 
smoothed  cumulative  distributions  of  Figure  26.  Therefore,  to 
the  extent  that  the  12  earthquakes  typify  the  earthquake  popu¬ 

lation,  and  the  11  surface  events  typify  their  population,  the 
curves  of  Eigure  27  may  be  taken  as  a  first  estimate  of  the 
probability  density  distributions  of  the  mean  cross-correlation 
of  P-coda  for  earthquakes  and  for  explosions.* 

The  surface  event  distribution  (h=0)  offers  no  surprises. 
However,  the  earthquake  distribution  is  unusual.  When  double¬ 
humped  distributions  of  this  kind  appear  in  nature,  it  is  often 
the  result  of  two  separate  phenomena  operating.  It  is  tempting, 
therefore,  to  break  the  distribution  of  Eigure  27  into  two 
roughly  symmetrical  distributions  as  shown  by  the  dashed  lines. 
It  is  further  tempting  to  ask  if  the  events  associated  with  the 
left-hand  part  of  the  earthquake  distribution  are  in  some  way 
distinguishable  from  those  associated  with  the  right-hand  part 
of  the  distribution.  We  find  that  the  left-hand  distribution 
Includes  all  the  earthquakes  deeper  than  the  Moho  while  the 

right-hand  distribution  represents  all  the  earthquakes  shallower 

than  the  Moho. 


*Here  it  is  necessary  to  caution  that  the  smallest  earthquake 
studied  had  a  magnitude  of  4.7. 
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This,  of  course,  is  very  attractive,  as  earthquakes  deeper 
than  the  Moho  can,  with  reasonable  confidence,  be  distinguished 
from  surface  events  by  other  means.  It  is  the  shallow  quakes 
which  constitute  the  problem,  and  in  Figure  27  their  distribution 
appears  separable  from  the  surface  event  distribution. 

On  the  assumption  that  the  deeper  earthquakes  can  be 
identified  by  other  means,  the  data  of  Figure  26  has  been 
replotted  as  Figure  28,  where  only  those  earthquakes  shallower 
than  40  kilometers  have  been  retained.  The  results  are  strik¬ 
ing.  Both  the  surface  events  and  the  shallow  earthquakes 
follow  Gaussian  distributions  as  closely  as  could  be  desired. 
Further,  the  separation  between  these  two  distributions  is  large. 
In  conventional  signal  processing  terms,  the  signal-to-noise 
ratio  is  16  dB. 

It  is  hardly  necessary  to  urge  caution  in  expressing  one's 
optimism  at  this  point.  The  foregoing  discussion  is  predicated 
on  too  many  if's  to  be  taken  seriously.  Consequently,  consider¬ 
able  effort  is  being  made  at  the  present  time  to  increase  the 
size  of  the  earthquake  population  manyfold.  Nevertheless,  one 
more  if  is  worth  observing;  mmely,  if  Figure  28  turns  out  to  be 
a  reasonable  approximation  to  the  real  world  distributions  of 
the  coda  correlation  for  shallow  earthquakes,  then  an  extremely 
powerful  discriminant  is  immediately  available. 
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FIGURE  26  -  FRACTION  OF  EVENTS  HAVING  MEAN  CODA  CORRELATION 


p  OR  GREATER 
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FIGURE  27  -  ESTIMATES  OF  THE  FIRST  PROBABILITY  DENSITY 
DISTRIBUTIONS  OF  MEAN  P-CODA  CORRELATION  OF 
SURFACE  EVENTS  AND  EARTHQUAKES 
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FIGURE  28  -  FRACTION  OF  EVENTS  HAVING  MEAN  CODA  CORRELATION 

Pc  OR  GREATER 


-53- 


GENERAL  ATRONICS  CORPORATION 


Report  1531-2033-7 


5  July  1966 


Section  IY  Attenuation  of  Seismic  Waves  in  the  Earth” s  Upper  Mantle 
General 

A  suggestion  by  Dr.  Eugene  Herrin  that  the  large  dis¬ 
crepancies  in  the  measured  spectra  from  LONGSHOT  might  result 
from  low  Q  regions  in  the  earth's  upper  mantle  prompted  a  brief 
study  to  discover  whether  it  is  possible  to  determine  the 
existence  of  such  layers  using  the  spectra  of  the  P  and  pP 
phases  of  deep  earthquakes.  This  section  describes  the  results 
obtained  to  date. 

The  absorption  strength  of  a  medium  can  be  measured  by  the 
fractional  energy  absorbed  per  cycle  =  Q.  Experiments  point 
toward  a  Q  constant  with  frequency  for  the  earth.  A  comparison 
of  the  frequency  content  of  two  phases  of  the  same  deep  event, 
pP  and  P,  makes  it  possible  to  find  the  average  value  of  Q  in 
the  upper  mantle  at  the  epicenter  region.  Three  events  from 
the  Fiji  Islands  region  were  analyzed.  These  events  yield 
measurements  of  Q  of  the  order  of  200,  indicating  that  the  atten¬ 
uation  is  one  order  of  magnitude  larger  than  expected  from  pre¬ 
vious  measurements.^ 

A  method  is  suggested  to  extend  the  study  to  treat  the  upper 
mantle  of  aseismic  regions. 

Description  of  the  Experiment 

The  attenuation  suffered  by  a  wave  traveling  through  a 
medium  is  characterized  by  the  dimensionless  factor  Q,  defined 
as 


=  fractional  energy  absorbed  per  cycle. 


2tt  _  AE 
Q  ~  E 

In  general  Q  is  a  function  of  frequency.  A  large  number  of 
observations  of  absorption  indicated  a  constant  energy  absorption 
per  cycle. ^  Since  no  linear  relation  between  stress  and  strain 
will  produce  a  Q  constant  with  frequency,  effect  other  than  simple 
absorption  must  contribute  to  the  empirical  results  obtained  by 
Howell.  This  fact  will  not  prevent  measurement  of  attenuation  on 
the  upper  mantle . 

For  a  given  range,  a  conslan"  factor  Q  implies  stronger 
attenuation  for  high  frequencies  than  low  frequencies.  The 
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attenuation  function  (energy)  for  a  monochromatic  wave  can  be 
considered  as 


A  (A,u,)  =  e-2a(u,)  * 

A  =  length  of  the  path 

a  (ou)  =  absorption  coefficient 

To  find  the  attenuation  function  that  corresponds  to  a  constant 
Q,  let 

A  =  X  (wave  length) 

AE  _  i  e-2aX  _  2tt 
E  Q 

2  2  2tt 

1  -  1  +  2aX  -  4a  X  + _ =  — 

Q 

for  aX  «1  :  a  =  —  = 

XQ  2cQ 

where  c  =  group  velocity 
~A 

A  (A,uu)  =  e”  =  attenuation  function  (energy) 

In  studying  attenuation  in  the  upper  mantle  the  phases  pP 
and  P  are  of  particular  importance.  The  path  of  pP  is  greater 
than  the  path  of  P.  This  additional  travel  is  composed  of  two 
crossings  of  the  upper  mantle  and  the  crust  in  the  region  of 
the  epicenter.  This  fact  allows  us  to  determine  the  attenuation 
characteristics  of  the  upper  mantle,  since  the  ratio  of  the 
energy  spectra  of  P  and  pP  is  determined  by  the  attenuation 
suffered  during  this  extra  path. 

Por  the  purpose  of  this  study,  the  following  model  will  be 
used  for  the  earth. 


P:  Pocus 

D:  Crust  thickness 
H:  Pocal  depth 
R:  Range  to  the  receiver 
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The  attenuation  functions  will  he  given  by  the  following 
approximate  expressions: 

,„)  =  e*V  (2'Ci) 


a) 


b) 


ApP  (R, 

=  length  of  path  in  region  i,  =  Q  value  in  region  i 

=  2(5^2)  +  —  +  (Terms  including  attenuation  in  the 
Gi  Q  lower  mantle  and  the  crossing  of 

the  upper  mantle  and  crust  at  the 
receiver  =  T) 


-s.  (2  (Ij2)  +  ££  +  I) 

ApP  (H,»)  =  e  0  Q  & 

Ap  (R,w)  =  e  c 


2D 


These  are  the  square  of  the  moduli  of  the  transfer  functions 
of  the  earth  as  a  filter,  for  both  cases,  and  therefore  relate 
the  energy  spectra,  at  the  receiver,  with  the  energy  spectrum 
at  the  focus. 

Por  the  following  two  cases  it  is  assumed  tta1;  the  source 
spectrum  for  P  and  pP  are  identical: 

a)  Spherically  symetric  radiation  pattern  (explosive  type 
of  source  mechanisms) 

b)  Plane  fault  as  the  source  mechanism,  P  and  pP  leaving 
the  source  approximately  in  the  same  direction  and  opposite  sense. 

Ep  (m,R)  _  Ap  M  _  +  |(2  (5j£)  +  ||) 

Epp(uj,R)  App(tu,R) 


where  Epand  E  p  are  the  energy  spectra.  This  means  that  we 
expect  to  fine!  more  attenuation  in  pP  because  it  crossed  the 
crust  and  the  upper  mantle  two  additional  times. 

Since  the  thickness  of  the  crust  is  small  (25-40  km),  the 
effect  of  the  attenuation  in  the  crust  (assuming  a  Q  of  200 
for  both  crust  and  upper  mantle)  is  very  small  compared  with  the 
effect  of  the  upper  mantle  crossing. 
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—  =  -—  =0.4  with  D  =  40  km 
Qc  200 

the  term  2(2=2)  =5.0  for  H-D  =  500  and  Q  =  200 
Q 

As  a  result,  should  the  average  Q  show  a  lov.  value,  this 
phenomenon  must  he  related  to  the  upper  mantle  and  not  the 
crust . 

Eor  an  earth  without  high  attenuation  in  the  upper  mantle, 
the  ratio  of  the  energy  spectra  should  have  the  following  expression: 


Ep  (in, R) 

Epp(n>,R) 


=  e 


ui/2h 

5 ' 


-u.  2Dx 


where  h  =  focal  depth 
^  =  average  Q  over  ray  path  =  2000 


2h 

t 

2D 

Qc 


1000 
2000 

^=  0.4 


=  .5 


In  this  case  the  attenuation  in  the  crust  could  not  be 
neglected;  however,  interpretation  of  the  results  would  be 
difficult  should  the  upper  mantle  have  an  intermediate  Q, 
(approximately  800).  It  was  the  purpose  of  this  investigation 
to  seek  only  existence  of  very  low  Q's  in  some  regions  of  the 
upper  mantle.  Hence,  we  can  side  step  this  problem. 

Experimental  Results 

Three  deep  earthquakes  in  the  Fiji  Island  region  were 
chosen  to  test  the  hypothesis  of  a  low  Q  layer  for  the  upper 
mantle  in  this  part  of  the  earth.  (See  Figure  29)  These  events 
were : 
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1)  E.Q . 

2)  E.Q. 

3)  E.Q . 


Fiji  Islands  Region  7  Oct  63 
25.6°S  179-9  E  d  -  550  *m. 

Fiji  Islands  Region  7  Dec  63 
22 . 1°S  179.4  W  d  =  -546  ion. 

Fiji  Islands  Region  26 April 64 
20 ,6°S  178. 0°W  d  =  490  km. 


13:14:24. 6z 
M  =  5.7 

04:07:52. Sz 
M  =  5.5 

14:52:07. 6z 
M  =  5.1 


The  energy  spectra  for  both  the  P  and  pP  phases  were  computed 
for  each  of  these  events,  and  the  ratio  of  the  spectra  of  the 
two  phases  was  then  compared  with  the  modeled  functions  for  a 
layer  with  high  Q  and  one  with  low  Q.  The  experimental  data 
together  with  the  theoretical  expressions  are  shown  in  Figures 
50  through  52  The  abscissa  is  frequency  in  cps.  The  ordinate  is 
the  logarithm  of  the  ratio  of  the  energy  spectra  normalized  to 
the  value  at  .5  cps.  Hence  all  curves  cross  the  abscissa  at 
0.5  cps.  From  these  curves  the  following  data  may  be  summarized. 
The  measured,  average  Q  of  the  upper  mantle  for  the  three  events 
studied  was  458,  218,  272.  As  can  be  seen  from  Figure  1,  the 
earthquake  leading  to  the  highest  measured  value  of  Q  is  the 
only  one  of  the  three  which  lies  on  the  Fiji  Ridge.  This  fact 
might  explain  the  difference  in  the  observed  values  of  attenuation. 
The  experimental  data  obtained  are  shown  in  more  detail  in  Figures 
50  through  52 .  These  values  show  a  tendency  to  have  a  smaller 
spread  with  respect  to  the  least  square  fit  as  we  increase  the 
record  length. 

The  least  square  fitted  curves  are  summarized  for  each 
event  in  Figures  50through  52 . 

Extension  of  the  Method  to  Analyze  Non-Seismic  Regions 

Some  seismic  stations  have  shown  a  greater  high  frequency 
content  in  their  records  than  others.  This  could  mean  that  the 
attenuation  characteristics  of  the  upper  mantle  are  regionally 
dependent.  To  analyze  this  hypothesis,  the  following  method  is 
suggested . 

Let  us  assume  that  we  are  interested  in  studying  the  atten¬ 
uation  characteristics  of  the  upper  mantle  of  a  certain  region 
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where  there  are  no  deep  earthquakes .  Let  us  call  that  area 
region  S.  We  must  first  study  some  other  region  which  possesses 
both  deep  earthquakes  and  a  seismic  station.  Let  us  call  this 
region  A.  Records  obtained  at  A  and  a  station  in  S  may  then  be 
compared,  taking  the  ratio  of  the  energy  spectra  of  the  P 
arrivals.  The  ratio  of  the  energy  spectra  of  the  P  phase  re¬ 
corded  at  A  and  S  enables  us  to  learn  about  the  attenuation 
characteristics  of  the  upper  mantle  in  the  region  S.  It  is 
important  to  note  that  in  following  this  procedure  care  must 
be  taken  to  include  any  frequency  dependent  radiation  pattern 
effects  which  will  result  from  geometry  of  the  measurement. 
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FIGURE  30  -  LEAST  SQUARES  FITS  TO  THE  EXPERIMENTAL  RATIO  OF  ENERGY  SPECTRA  (SUMMARY) 


REGION:  FIJI  ISLANDS 

SEIS.  NO.:  5296 

DATE:  7  DEC  63 

ORIGIN  TIME:  4:7:52.8 

DEPTH:  546  KM 
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FIGURE  31  -  LEAST  SQUARES  FITS  TO  THE  EXPERIMENTAL  RATIO  OF  ENERGY  SPECTRA  (SUMMARY) 


REGION:  FIJI  ISLANDS 

SEIS .  NO.:  2106 

DATE:  26  APRIL  64 

ORIGIN  TIME:  14:52  =  07-6 

DEPTH:  490  KM 


(UamtAL  ATKOJUUS  CORPORATION 


(SdO  5 


2 

X 
X 
Eh  Eh 

O 

CO  X 

H 

X  < 

CO 

X  s 

CO 

o 

X 

X  X 

X 

O  X  X 

X 

H  X  X 

X  H 

CO  EH  X 

X 

CO  X 

11 

X  X 

Q 

X  H  * 

X  Q 

X  o 

Eh  \ 

\  X  o  o 

O  X 

\  x  o  o 

X 

1  CM  CM 

X  X 

\  Q 

X  o 

\  X  11  II 

X  X 

1  Eh 

\ 

EXPEC 
FOR  Q 
Q 

o 

CM 


l n 


o 

x 

00 


US 

Eh 

U 

2 

X 

X 


o  o 
o-o 

H  O 

i — i 


x  x 

Oh  Oh 


in  o 

CM  O 


>H 

X 

< 

s 

X 

CO 


< 

X 

Eh 

o 

W 

x 

CO 


>H 

o 

X 

x 

x 

x 

o 

o 

H 

Eh 

< 

IX 


< 

Eh 

2 

X 

X 


H 

X 

x 

Oh 

X 

w 


w 

X 

Eh 


O 

Eh 


CO 

Eh 

H 

X 


CO 


< 

X 

O' 

CO 


Eh 

CO 

< 

X 

X 


I 


Q  HO 
X  i — I 

o 

o  ••  •• 
X  XX 
X 

•  • 


CM 

on 

X 

X 

x 

o 

H 

X 


o  o 

CM  «H 

0  01  aaznVWHON  OliVH  Vh£o3<3S  XDH3N3)  dot  ot 

-63- 


o 

on 


GENERAL  A2R0NICS  CORPORATION 


(sjo  s*o  01  aaznvwHON  oiivh  vnioaas  xonaNa)  dot  ox 


64 


FIGURE  33  -  ENERGY  SPECTRA  RATIO  -  EXPERIMENTAL  DATA  AND 

LEAST  SQUARES  FIT 


REGION:  .  FIJI  ISLANDS 

SEIS .  NO. :  53dl 

DATE:  7  OCT  63 

ORIGIN  TIME:  13:l4:jl).6 
DEPTH:  550  KM 
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FIGURE  31*  -  ENERGY  SPECTRA  RATIO  -  EXPERIMENTAL  DATA  AND  LEAST  SQUARES  FIT 


REGION:  FIJI  ISLANDS 

SEIS .  NO.:  5301 

DATE:  7  OCT  63 

ORIGIN  TIME:  13:14:24. b 
DEPTH:  550  KM 
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FIGURE  35  -  ENERGY  SPECTRA  RATIO  -  EXPERIMENTAL  DATA  AND  LEAST  SQUARES  FIT 
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FIGURE  36  -  ENERGY  SPECTRA  RATIO  -  -  EXPERIMENTAL  DATA  AND  LEAST  SQUARES  FIT 


REGION:  FIJI  ISLANDS 

SEIS .  NO.:  5296 

DATE:  7  DEC  63 

ORIGIN  TIME:  A : 7 : 52 . 8 

DEPTH:  5^6  KM 
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FIGURE  37  -  ENERGY  SPECTRA  RATIO  -  EXPERIMENTAT  LATA  AND  LEAST  SQUARES  FIT 
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FIGURE  38  -  ENERGY  SPECTRA  RATIO  -  EXPERIMENTAL  DATA  AND  LEAST  SQUARES  FIT 
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FIGURE  39  -  ENERGY  SPECTRA  RATIO  -  EXPERIMENTAL  DATA  AND  LEAST  SQUARES  FIT 
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FIGURE  40  -  ENERGY  SPECTRA  RATIO  -  EXPERIMENTAL  DATA  AND  LEAST  SQUARES  FIT 


REGION:  FIJI  ISLANDS 

SEIS .  NO.:  2106 

DATE:  2*1  APR  6*1 

ORIGIN  TIME:  1*1:52  =  07-6 

DEPTH:  *190  km 
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FIGURE  *11  -  ENERGY  SPECTRA  RATIO  -  EXPERIMENTAL  DATA  AND  LEAST  SQUARES  FIT 
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